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SUMMARY 


Higher order compact algorithms are developed for the numerical simulation of wave propagation by 
using the concept of a discrete dispersion relation. The dispersion relation is the unique imprint of any 
linear operator in space-time. The discrete dispersion relation is derived from the continuous dispersion 
relation by examining the process by which locally plane waves propagate through a chosen grid. The 
exponential structure of the discrete dispersion relation suggests an efficient splitting of convective and 
diffusive terms for dissipative waves. Fourth- and eighth-order convection schemes are examined that 
involve only three or five spatial grid points. These algorithms are subject to the same restrictions that 
govern the use of dispersion relations in the construction of asymptotic expansions to nonlinear evolu- 
tion equations. A new eighth-order scheme is developed that is exact for Courant numbers of 1, 2, 3, 
and 4. Examples are given of a pulse and step wave with a small amount of physical diffusion. 


INTRODUCTION 


Numerical methods for computing dynamic phenomena, especially wave propagation, are becoming 
a major focus of modem research because of the pervasive nature of unsteady flows in all areas of fluid 
mechanics and the increasingly wide availability of supercomputer resources. The propagation of waves 
over long distances requires precise phase accuracy. It is an unfortunate property of all finite-difference 
schemes that phase resolution must be compromised over some portion of the available spectrum. In 
order to correctly capture wave phenomena, the spatial and temporal order of the discretized equations 
must be accurately matched to a relatively high order. Otherwise, extremely small time steps must be 
taken to maintain phase resolution, which can result in the buildup of roundoff error and the inefficient 
use of computational facilities. 

In general, currently available algorithms for wave propagation may be placed in two major cate- 
gories: Crank-Nicolson (C-N) and characteristic type. The former approach was originally developed in 
reference 1 for parabolic equations and extended by others to hyperbolic problems, with mixed success. 
(The use of the term Crank-Nicolson for implicit finite-difference methods applied to nonparabolic 
equations follows the nomenclature adopted by Mitchell in ref. 2.) The characteristic-based method 
(apparently first described in ref. 3) was developed especially for wave-like phenomena. 

The C-N approach is a very successful strategy to attain high accuracy in a given computational 
molecule by using implicit finite differences. This method has withstood the test of time and is clearly 
the method of choice for highly dissipative systems. The early promise of using central differences for 
convective systems was short-lived (ref. 4). It is now apparent that certain low order central-difference 



algorithms are not useful for resolving convection-dominated flows. This class of problems is the essen- 
tial feature of many practical wave problems where the ratio of the convective to diffusive terms can be 
as large as 10 8 . 

Characteristic-based differencing was first presented in the explicit first-order scheme of Courant, 
Issacson, and Rees in 1953 (ref. 3). This is equivalent to tracking the retrograde characteristic with a 
linear interpolation between grid points. Since then, a wide variety of these upwind algorithms have 
been developed. A major breakthrough in the analysis of upwind methods was the development of 
so-called Godunov schemes. These algorithms model the equations of convective gas dynamics with a 
special treatment of the discontinuity surfaces by solving a sequence of Riemann problems (called the 
“breakdown formulas” in ref. 5) and imposing a special monotonic criterion on the wave. This concept 
was exploited in a large number of modem unsteady algorithms introduced in the early 1980s. A review 
article concerning these methods as applied to inviscid gas dynamics was recently written by Roe 
(ref. 6). These algorithms seem to be the most evolved versions of the simple upwind schemes. In prac- 
tical applications they can be very complicated for waves traveling in multiple directions and for the 
inclusion of viscous effects. 

This paper describes a different approach for generating finite-difference approximations to 
convection-diffusion problems. The response of the operator to sinusoidal wave trains is considered, 
rather than attempting to deal with local characteristics or Taylor-series expansions. The goal of the 
method is to develop the best possible approximation to the space-time differential operator, with the 
constraint that the system of discretized algebraic equations maintain a specified minimal bandwidth. 
The measure of the “best possible approximation” is that the discrete transfer function (also called the 
amplification factor) that advances the solution at each time step has minimal dissipative and dispersive 
error. The key to generating the optimal finite-difference equation is to carefully examine the dispersion 
relation, which is an intrinsic property of every differential operator. In particular, the linearized theory 
of dispersive waves has been developed to a very high order (refs. 7 and 8). With this theory, algorithms 
are developed using the analogy between asymptotic dispersive-wave propagation and local plane-wave 
propagation through a discrete mesh. 

In the following sections, the discrete dispersion relation (DDR) is defined and applied to some well 
known algorithms. The DDR also suggests an efficient time- splitting to separate convection and diffu- 
sion in one dimension, or even to split multidimensional problems The important role of the data band- 
width is illustrated by the propagation of pulses with small diffusion. The propagation of a step function 
is compared with the monotone-preserving scheme of Leonard (ref. 9), where it is shown that the 
“wiggles” can be effectively suppressed with small values of physical diffusion in the differential 
operator. 

The method described here is contrasted with another class of algorithms called compact schemes 
(or Hermitian methods in ref. 10). These methods are based on the fact that a differential equation is 
satisfied at each of a number of mesh points in the chosen computational molecule, thus relating the 
dependent variable and its derivatives at these spatial mesh points. Such methods have been used with 
considerable success for convection-diffusion problems. The new method attempts to generate similar 
compact schemes by exploiting the differential operator in space-time rather than by treating only the 
spatial response. 
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A fourth-order scheme is derived that is exact for values of the Courant number (Cn) equal to 1 and 
2, and an eighth-order scheme that is exact for Cn equals 1, 2, 3, and 4. These should be very efficient 
for many wave problems, since no decision regarding the local velocity direction need be made. The fact 
that only three (or five) mesh points are needed for the fourth- (or eighth-) order schemes should make 
them very useful. In theory, although the method is neutrally stable for all values of Cn, there is a prac- 
tical limit to this parameter because of severe wave distortion at short wavelengths, due to increasingly 
large dispersive errors. 

Recently, a series of higher order algorithms of the C-N type have been developed by Noye and 
co-workers (refs. 1 1 and 12). They used the method of modified partial differential equations (ref. 11) 
and Pade approximations (ref. 12), following the ideas presented by Young in the appendix of refer- 
ence 13. The developed schemes depend on a variety of parameters, and Noye has also derived an 
“optimal” fourth-order scheme for convection that is the same as the one developed here using the DDR. 
In reference 12 a seventh-order scheme is derived, based on the Pade method. 

In the next section, the concept of the DDR associated with the convection-diffusion operator is 
derived and compared with other approaches. Although Fourier methods have been used extensively in 
the past to examine the stability and dispersion of numerical approximations, they have not been used as 
the primary tool to derive computational molecules. Later sections will compare results for a number of 
model problems regarding the convection and diffusion of pulses and steps. Finally, a description of 
possible future work will be presented. 


THE DISCRETE DISPERSION RELATION 


The basis for developing high Fidelity Finite-difference equations is to examine the manner in which 
the locally linearized equations propagate plane waves. The equation satisFied by the wave number 
determines a dispersion relation at the chosen time step. This dispersion relation may vary from point to 
point, and at each time step, for systems with time-dependent coefficients. The associated difference 
equation is found from a local Maclaurin series which relates the exact dispersion relation to the 
assumed approximation at the desired order. 

First, consider the prototype one-dimensional convection equation 

du ,,du _ 

■5T + V 5x = 0 (1) 

A plane-wave solution of the form 

u(x,t) = Ae i(ft,t+kx) (2) 
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is substituted in equation (1) to obtain the simple dispersion relation 


0) + Vk = 0 


( 3 ) 


If x is a discrete time increment, the ratio of the exact solution at time nx and (n+l)x at the point x is 

n+1 i[co(n+l)x+kx] 
u e icot 

u n e i(an+kx) 


(4) 


which does not depend on x and is independent of n. (This ratio is homogeneous in space and station- 
ary in time.) If the spatial dimension is discretized with a uniform step size h, the dispersion relation in 
equation (3) is used to relate the time step to the spatial scale, using the Courant number (Cn = Vx/h) as 
the scaling parameter. The final result is 


,n+l 


= P = e 


-iCn kh 


(5) 


Equation (5) relates the transfer of information to the next time step at a given point in terms of the 
spatial scales. The symbol p represents the exact amplification factor on a finite grid superimposed on 
the space of independent variables. This spatial-temporal coupling is embedded in the dispersion rela- 
tion. In this extremely simple example, the relation is purely linear with the Cn the constant of 
proportionality. 

Commonly used backward and central finite-difference schemes on a uniform mesh can also be 
expressed as amplification factors. These are approximate transfer functions that relate to equation (5) 
only in terms of matching its Maclaurin series in kh to a given order. Amplification factors for two well 
known schemes are 


„n+l 

- — = (1-Cn) + Cne _lkh 


(6) 


^ = (1 - Cn 2 ) + icn(l + Cn)e -ikh - ^Cn(l - Cn)e ikh 


The first of equation (6) is the Courant Issacson and Rees method of reference 3, while the latter is the 
Lax-Wendroff scheme. The mesh can support all wavelengths up to 2h, so kh must vary between 0 
and 7t. 
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The most general three-spatial-point/two-time-level stencil is 



( 7 ) 


The dependent variable at the discrete time n + 1 and centered at discrete point j is expressed in 
terms of its nearest neighbors at times n and n + 1. The fundamental question is not how to determine 
the coefficients from a predetermined spatial/temporal finite-difference formula, but how to use the 
available constants to best approximate the underlying partial-differential operator, in this case, equa- 
tion (1). A large number of implicit/explicit and/or biased/unbiased schemes can be constructed using 
the available constants. 

The computational molecule described by equation (7) has a local influence property that is an 
analog to asymptotic wave theory. Thus, in reference 14 a wavetrain is defined as “a system of almost 
sinusoidal propagating waves with a recognizable dominant local frequency and wave number.” This 
wave train may be propagating in a nonhomogeneous medium, but the local wave properties depend on 
the dispersion relation, which is defined in terms of the locally constant coefficients of the operator. The 
same statement may be made with respect to the six-point computational molecule in equation (7) where 
a “test wave” of local frequency co and wave number k propagates through the six-point computational 
molecule. Determining the coefficients now becomes a matter of mapping the dispersion relation — the 
unique imprint of the differential operator — onto the grid, and using this relation to generate as accurate 
an amplification factor as possible. 

Let a locally plane wave of the form e icot+lkx propagate through the mesh. The amplitudes are 
related by equation (7): 

e i(o(n+1)x (a 0 + aj t~ iVh + a 2 e ikh ) = e iajnT (b 0 + b, e" ikh + b 2 e** 1 ) (8) 


The amplification factor (3' is the ratio of the solution at subsequent time steps: 


Jot _ bp + b i e ^ + b 2 e lkh 
13 = 6 "ao + a 1 e- ikh + a 2 e ikh 


(9) 


The exact ratio is given by (3 in equation (5). It lies on the unit circle in the complex plane for all 
values of Cn and kh, since the underlying operator is dissipation-free. The approximate operator in 
equation (9) can be made free of dissipation if the numerator and denominator are complex conjugates to 
one another. This condition is strictly enforced if ao = bo, t ai = b 2 5 and &2 = bi. The expression for (3 
now contains three constants, of which only two are independent. The only error that remains is that due 
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to phase dispersion. Taking account of the fact that the numerator and denominator are complex conju- 
gates of the form (A + iB)/(A - iB), a simple form for the argument of Num[P'] is 


tan[arg(Num[P'])] = -5- 


( 10 ) 


From equation(4), the exact value of the argument is 


tan 




( 11 ) 


The factor 1/2 is used because only half the phase shift is attributed to the numerator of the discrete 
amplification factor. 

The available constants can be used to match these equations up to order (kh)^. The Maclaurin series 
is derived in the appendix where the coefficients are specified as a function of Cn. The final formula is 

^(Cn - l)(Cn - 2)uj*_4 - i(Cn - 2)(Cn + 2)uJ +1 + ^ (Cn + l)(Cn + 2)uJ 4 ? 

= -j^(Cn + l)(Cn + 2)uj_! - ^(Cn - 2)(Cn + 2)uJ + j^(Cn - l)(Cn - 2)uJ +1 (12) 


This is a fourth-order accurate approximation to equation (1) where the order is defined as matching 
the Maclaurin series to the indicated number of terms. This is the formula derived by Noye in refer- 
ence 11, using the modified equation approach. His method used weighed finite-difference formulas in 
both space and time to obtain higher order schemes. His “optimum” scheme is the one shown in equa- 
tion (12), which is the unique result of approximating the DDR with a three-spatial-point/two-time-level 
scheme. 

Exactly the same procedure can be used to generate higher order dissipation-free formulas. If one is 
willing to accept wider bandwidths, arbitrarily high orders can be generated. Such formulas may be 
considered high fidelity versions of the C-N method. For example, let a locally plane wave propagate 
through the five-spatial-point/two-time-level stencil. The resulting eighth-order scheme (also derived in 
the appendix) is 
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(Cn - l)(Cn - 2)(Cn - 3)(Cn - 4)u] 1 _ + 2 1 - ^ (Cn - 4)(Cn - 3)(Cn - 2)(Cn + 4)u] , _ + 1 1 
+ 2^o (Cn - 4)(Cn - 3)(Cn + 3)(Cn + 4)uJ +1 - ^ (Cn - 4)(Cn + 2)(Cn + 3)(Cn + 4)uJ + + 1 1 

+ mo' (Cn + 1)(Cn + 2)(Cn + 3)(Cn + 4 )u j ++ 2 

= (Cn + l)(Cn + 2)(Cn + 3)(Cn + 4)uJ_ 2 - 4^1 (Cn - 4)(Cn + 2)(Cn + 3)(Cn + 4)uJ_i 
+ 2J0 (Cn - 4)(Cn - 3)(Cn + 3)(Cn + 4)uJ - 4^' (Cn - 4)(Cn - 3)(Cn - 2)(Cn + 4)uJ +1 

+ T^0 (Cn " 4)(Cn “ 3)(Cn " 2)(Cn " 1)U J+ 2 (13) 

A direct extension of the second-order C-N scheme from parabolic to convective equations can be 
obtained in a simple approximation to equation (7) with bo = ao = 1 and bi = -b2 = -ai = a 2 . This 
formula is well known and is derived in another manner by Mitchell (ref. 2, p. 167). 

The local order of accuracy of a given finite-difference approximation is defined as the power of the 
first nonvanishing term in a Taylor series expansion of the difference between the exact and discrete 
operators. Since the underlying dispersion relation itself represents the dominant term in an asymptotic 
expansion for nonuniform or nonlinear wave propagation, the formulas given above may not have the 
same local order of accuracy with respect to the differential operator. In fact, it can be shown that the 
schemes are of fourth- or eighth-order accuracy only if the underlying equation has constant coefficients. 
This will be illustrated with fourth-order method in equation (12). 

Expand each of the terms in equation (7) in a Taylor series about the point u" with spatial and tem- 
poral increments h and x. Using the coefficients in equation (12) and after some complicated algebra, 
the final form is 

u t + Vu x = ^ - ^x(u tt + Vu xt ) - ^-t 2 (4u ut + 6Vu xtt + 2V 2 u xxt ) 

-34 ' * 3 (utttt + 2Vu xtu + V 2 u xxtt ) - ^h 2 (u xxt + Vu xxx ) 

h 2 x(u xxtt + Vu xxxt ) + 0(x 4 , h 4 ,x 2 h 2 ) (14) 

where 

A(u) = aoUj* +1 + aiUj 1 ^ 1 + a 2 uj l 4 !| 1 - bou" - bjUjLi - b 2 Uj + j 
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If the convection velocity V is constant all of the error terms above vanish and the discrete operator 
is a fourth-order accurate space-time approximation to both the dispersion relation and the differential 
equation. If V is spatially variable, only the time derivatives can be factored and the local order of 
accuracy drops to second order. Finally, if V is variable in space and time, the local order of accuracy 
drops to first order. This behavior is consistent with the concept of a dispersion relation as an exact rep- 
resentation for solutions of partial differential equations of simple form (constant coefficients), but only 
as the leading term of an asymptotic series for equations of more complex structure. In the past, success- 
ful use of the continuous dispersion relation to explain major physical features of wave propagation 
indicates that the use of the DDR may be useful in selected simulations. For example, the propagation of 
acoustic or electromagnetic waves usually involves systems with constant coefficients, and the analysis 
of linear stability in fluids involves wave motion governed by systems with spatially variable coeffi- 
cients. In such cases, the resolution of the numerical solution may depend more on the stability proper- 
ties (that is the amplification factor) of the discrete operator than on its rate of convergence as reflected 
by the local order of accuracy. 


Now consider the effect of dissipation using the second-model differential operator. 


5u 

di 


+ V 


3u 

9x 



(15) 


Solutions in form of wave trains u = A eKwt+kx) yields the dispersion relation for waves in a diffu- 
sive medium 


id) + iVk + vk 2 = 0 


(16) 


where V and v are real constants. This equation relates the real wave number k to the radian frequency 
CO in terms of the real parameters V and v. A nondimensional form of equation (15) in terms of the 
mesh parameters is 


cot * -Cn kh + iDn(kh) 2 


(17) 


where the Diffusion number, Dn, is defined as vx/h 2 . Following equation (5), the exact dispersion rela- 
tion is 


P _ e ion _ e -iCnkh-Dn(kh) 2 



(18) 


which can be split into the product of convective and diffusive terms, as shown. Since the dispersion 
relation is derived from a differential operator with locally constant coefficients (following its use in 
wave theory), it will always support waves of exponential form. Furthermore, the number of terms in the 
exponent — the second exponential in equation (18) — will always sum to the number of space-derivative 


8 


terms in the operator. Thus, the dispersion relation can always be broken into a product of transfer func- 
tions; in this case, one amplification factor for convection and one for diffusion. DDRs for convection 
have already been considered. The discrete transfer function for dissipation is not so critical for wave 
processes. It is the nature of dissipative processes always to evolve toward the low end of the wave- 
number spectrum. If the usual C-N algorithm is used as the dissipative discrete transfer function, a 
complete algorithm for convection and diffusion is obtained. 

The convective DDR generates an intermediate solution that is further processed by the diffusive 
DDR to obtain the complete solution at the next time level. Boundary conditions for the intermediate 
solution are of the convective type, and are almost always obvious for waves in infinite media. The same 
arguments can be used to generate a sequence of transfer functions to treat other spatial directions, or 
even terms without derivatives (such as arise in the Helmholtz equation). The use of the dispersion rela- 
tion is the unifying principle in generating approximate solutions to complex systems as a sequence of 
simpler operators. 


EXAMPLES OF THE DISCRETE DISPERSION RELATION 


The DDR has been shown to be a useful tool for generating approximations to differential operators, 
and for separating convective and diffusive effects at the algorithmic level. In this section some repre- 
sentative algorithms will be presented and contrasted with the fourth- and eighth-order convection meth- 
ods derived in the previous section. It is found that most algorithms do not perform very well for Cn 
greater than one. Even if they are stable, they exhibit very large phase dispersions that make them less 
useful for unsteady flow analysis. 

Results are presented for a range of Cn greater than unity. Many commonly used explicit methods 
are unstable for Cn > 1, but for comparative purposes three available algorithms are considered. Fig- 
ure 1 shows the magnitude and phase characteristics of the amplification factor for a first-order upwind 
scheme (ref. 3 ), a second-order Lax-Wendroff scheme, and a third-order backward scheme due to 
Leonard (ref. 15). Curves are compared with the theoretical dispersion relation from equation (5) for Cn 
from 0.5 to 3.0 in steps of 0.5. Notwithstanding the fact that that all methods are unstable for Cn > 1, 
attention is focused on the phase. The first- and third-order methods have the interesting property that 
the phase response is correct for Cn = 0.5. All methods are exact for Cn = 1.0 in both amplitude and 
phase, and the third-order method is exact for Cn = 2. The second-order method has the distinct advan- 
tage that no decision regarding the sign of the velocity is needed, and only a three-point stencil is 
needed. 

Figure 2(a) shows the same quantities for the fourth- and eighth-order methods. The amplitude is 
unity for all values of Cn and the fourth-order method is exact for Cn = 1 and 2. Figure 2(b) represents 
the dispersion characteristics for the eighth-order method over a larger range of Cn. It is exact for 
Cn = 1, 2, 3, and 4. 

This increased accuracy does not depend on an ever-widening computational molecule. The fourth- 
order method uses the same three-point stencil as the Lax-Wendroff method and the eighth-order method 
uses the same five-point mesh as Leonard’s third-order algorithm. The increased accuracy is achieved by 
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examining the differential operator (in terms of its dispersion relation), rather than the space and time 
derivatives separately. 

A measure of the phase accuracy is the quantity 0d, the departure wave number as defined by the 
value of 0 where the DDR and the exact dispersion relation first diverge. Since the numerical phase at 
kh = jt must always be zero or an integer (see fig. 2), 0d must oscillate as Cn increases, with a maxi- 
mum value of it when the algorithm is exact. For the eighth-order scheme, 0d,min is about 0.771 and 
decreases rapidly for Cn > 4. In contrast, the phase departure for the fourth-order scheme has a mini- 
mum of about 0.57C and decreases rapidly for Cn > 2. The important factor of phase fidelity is not 
always appreciated, and many common algorithms depend on numerical dissipation (“artificial viscos- 
ity”) to avoid oscillations. Such algorithms that claim no limits on the Cn based on stability considera- 
tions must be carefully examined to ensure that the numerical dissipation not only suppresses spurious 
oscillations, but does not damp out important high-wave-number physical phenomena. 

Another way of comparing phase error is shown in figure 3. Curves of phase error, defined as the 
ratio of the approximate to the exact phase shift, are drawn on a polar diagram with the phase kh as 
polar angle. The unit half-circle defines the locus of perfect phase resolution. A separate plot is shown 
for three values of Cn: 0.2, 0.8, and 4.0. Each plot shows the phase resolution for three dissipation-free 
algorithms: the second-order method reported in Mitchell (ref. 2), and the fourth- and eighth-order 
methods described above. The higher order algorithms follow the unit half-circle for a much greater 
length of arc, indicating higher phase fidelity. (For reference, the ray corresponding to kh = ji/ 2 corre- 
sponds to a disturbance wavelength of 4h.) The plot for Cn = 4 is one of the special cases where the 
eighth-order method is an exact analog of the differential operator (when V = constant). The phase 
resolution for the other methods is very poor. 

The effect of phase fidelity on the numerical simulation of convection operators, and the use of 
physical diffusion to supresses oscillations, will be examined in the following section where the resolu- 
tion of pulse and step-function wave motion will be examined. 


COMPUTATIONAL EXAMPLES 


The following examples are based on numerical solutions to equation (15) using the finite-difference 
approximations given by equations (12) and (13). Both examples are meant to illustrate the wave- 
capturing capabilities of the method and the importance of phase accuracy. The examples deal with 
propagation in an infinite medium and, where necessary, periodic boundary conditions are used. 

A feature of numerical wave propagation with Cn greater than one is that the system of algebraic 
equations may become ill-conditioned or the coefficient matrix may even be singular. In most cases 
periodic boundary conditions assure diagonal dominance, and this was found to be true for the systems 
generated by equations (12) and (13). 
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Example 1. Pulse Propagation in an Infinite Medium 


The first example will deal with the propagation of a dissipative pulse. This problem illustrates the 
critical role that the amplification factor plays as a filter through which the solution is processed at each 
time step. It will be demonstrated that the bandwidth of the data is the main parameter governing the 
accuracy of the simulation. A test pulse with the largest possible bandwidth is a Dirac delta function, an 
unattainable ideal in numerical simulations. The nonlimiting form of the Dirac function (ref. 16) is an 
exponential function of the form 


u(x,t) = 



(19) 


where n = t^/4(t - t r ), X = x - x r - t/tc and x,t are points in space and time. The Dirac delta function 
would correspond to the limit n — » 

The available physical parameters are x r the initial location of the pulse; t r the time at which the 
pulse is singular (t r is negative); tc the convective time scale; and t<j the diffusive time scale. (With 
respect to equation (15) and a length scale L, tc and td are L/V and L 2 /v, respectively.) A wave with a 
large ratio of td/tc would represent a convection-dominated flow. (The ratio td/tc is the Reynolds 
number.) Very small values of the diffusive time scale would represent a highly dissipative wave, and 
large values of tc would correspond to a sharp initial pulse with a rich frequency content. 

Equation (19) has the agreeable property that its Fourier transform with respect to X retains a simi- 
lar functional form. A numerical solution to equation (15) is generated using u(x,0) as initial conditions. 
The two new parameters that arise from the discretization are Cn and Dn, as described previously. 
Results will be shown for both the fourth- and eighth-order convection algorithms. 

A comparison for a relatively narrow banded pulse is shown in figure 4. The pulse has convected 
without a change in form, as shown in figure 4(a). Figure 4(b) presents the magnitude of the finite 
Fourier transform on an abscissa that is matched to the quantity kh/7t. The bandwidth of the pulse is 
approximately 0.6 tc. Figures 4(c) and 4(d) show the phase characteristics of the DDR as compared to the 
exact phase (dashed). Since the phase departure wave number is located well above the pulse bandwidth, 
as shown in figure 4(b), both the fourth- and eighth-order algorithms give similar predictions. The two 
computed waveforms are indistinguishable from the exact theoretical values, all three of which are 
shown superimposed as pulse B in figure 4(a). 

Figure 5 shows a similar computation for a large increase in the Cn to 4.78. This value corresponds 
to almost 2.5 revolutions of the DDR about the origin in the complex plane and is a relatively severe test 
of the phase resolution. The predicted final pulse waveform for the fourth- and eighth-order schemes are 
shown in figure 5(a). The fourth-order algorithm exhibits a strong lagging oscillation while the eighth- 
order method still correctly tracks the pulse. Since the differential operator, as well as the DDR, is 
dissipation- free at each time step, the magnitude of the amplitude spectra as shown in figure 5(b) are 
matched for both schemes, even though the phase of the fourth-order method does not accurately span 
the spectrum of the disturbance. What happens is that the real and imaginary components of the fourth- 
order spectrum have been severely redistributed above 0(j. This is shown in figure 6, which compares 
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the real and imaginary parts of the spectrum with the exact values. The harmonic content of the compu- 
tation (solid line) leads the the exact value (dashed line) for those wave numbers greater than 0d- Thus, 
the time-domain pulse shows a severe phase redistribution which results in a train of lagging “wiggles.” 

This result illustrates the strong problem-sensitive nature of the finite-difference formulation. In 
some applications, it may be possible to capture the correct partition of energy (represented by the spec- 
tral magnitude) in wave number space, even though significant phase errors occur. Other applications, 
such as acoustic wave propagation, require accurate phase tracking throughout the computation. 

The next result shows the case where the spectral bandwidth is increased by choosing a sharper 
initial pulse — the value of n in equation (19) goes from 5,000 to 10,000. A series of calculations that 
include dissipation are shown in figure 7. These solutions are obtained in the manner described in the 
previous section, where the convective and diffusive terms are separately modeled with individual dis- 
crete dispersion relations, each having its own finite-difference formulation. The solution is obtained as 
a two-step process. The exact solution is shown as a dashed line and is almost coincident with the com- 
puted solution that is shown as connected points. Finally, figure 8 shows the DDR and magnitude for the 
intermediate value of Dn. The slight dissipation affects the amplitude of the pulse and is reflected in a 
spectrum that evolves to lower wave numbers. The damping is not quite sufficient to supress the lagging 
oscillations shown in figure 7 (middle panel). 

In this simple first example, the spectrum was fixed by the initial data. The next case will show the 
much more severe test when a step function must propagate through the mesh. 


Example 2. Propagation of a Unit Step Function 

The previous section has indicated how the bandwidth and dissipative properties of the DDR can be 
tailored to correctly resolve a band-limited wave. Such flexibility is not possible in the general case, and 
an extreme example is the propagation of a finite discontinuity. A simple step function contains impor- 
tant phase and amplitude information at all wavelengths. Unlike the Fourier transform of the Heaviside 
Step Function in an infinite domain, the finite Fourier series is more complicated. 

oo 

f(x) = x i 2 V sin(n7txo)cos[2n7c(x-x 0 /2)] ^0' 

W — x 0 2 ^ n 

n=l 


which represents 


f(x) = l 0<x<XQ 

f(x) = 0 xq<x<1 

f(x + l) = f(x) 
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The harmonic amplitudes decay as 1/n, which is not fast enough to neglect higher harmonics. In 
addition, the power spectrum has a regular progression of zeros at the nodes of the sine function, which 
results in a highly regular series of hills and valleys. This feature of the spectrum, which is typical of 
propagating discontinuities, must be very well resolved by the DDR at all wavelengths in order to have a 
high fidelity simulation. This harmonic richness makes the step function an ideal model with which to 
study the performance of numerical algorithms. Leonard (ref. 9) has recently published a detailed study 
of the propagation of a step function as governed by the model convection equation (1). In this section, 
the convection of a step will be compared with Leonard’s results. 

The relative dispersion of second-, fourth-, and eighth-order schemes are shown in figure 9 for a Cn 
of 0.5. Directly above these curves is the magnitude spectrum for a step function from equation (20) on 
the unit interval. (Note that a logarithmic scale is now used for the ordinate.) In order to propagate the 
step with high fidelity, the phase shift at each wave number must increase linearly, as shown in the 
dashed line. Finite difference formulas are incapable of such fidelity, since the phase shift at kh = n is 
constrained to be an integer (in this case, zero). For this reason, the step is always distorted and wiggles 
always appear. The amplitude and dominant frequency of the oscillation, however, can be controlled 
with higher order schemes. Figure 10 shows a step function starting from xo = 0.22 after 45h intervals 
for Cn = 0.5. As the order of the underlying algorithm increases, the amplitude of the wiggles decrease 
and their dominant frequency increases. The traditional method of controlling this numerical artifact is 
to use a quantity called the “artificial viscosity.” In upwind schemes, even with forms of higher order, 
this artificial viscosity is built into the algorithm. This yields smoother solutions, but may result in unac- 
ceptable smoothing in some cases. 

In the current example, physical diffusion as it appears in equation (16) is used to induce just enough 
damping to suppress the oscillations. This approach has the advantage that a quantitative lower limit of 
the physical viscosity may be found and used to assess the practicality of a particular algorithm. The 
effect of a small amount of physical viscosity is shown in figure 1 1 for the fourth-order scheme and in 
figure 12 for the eighth-order scheme. 

These figures show that higher order schemes require smaller diffusion numbers to suppress oscilla- 
tions. This increases the possibility for high definition solutions within available mesh and time-step 
constraints. The trade-off among the equation bandwidth, algorithm complexity, and available storage 
must depend on the particular physical problem. 

The effect of diffusion on sharpening up the pulses is shown in the amplitude spectra in figure 1 3 for 
the eighth-order scheme. The diffusive factor in the DDR smooths out the high wave number as indi- 
cated in the figure. The predominant wavelength of the oscillation for the undamped wave is evident 
from the mismatch between the hills and valleys from the exact spectrum (dotted line) and the computed 
spectrum (dashed line). The most severe mismatch is at the departure wave number 0d > as shown in fig- 
ure 9. If sufficient damping is included to just cancel out this critical mismatch, a high resolution simula- 
tion can be achieved. Too much damping will detune the simulation, resulting in an over-smoothed dis- 
continuity, as is achieved with first-order upwinding. 

A summary chart of the step function results is shown in figure 14. This figure depicts the error as 
abscissa and the the variation of the error (wiggles) as ordinate. The region to the lower left is most 
desirable, the region to the lower right is the smoothed upwind type of solution, and the solution along 
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the diagonal is the “wiggly” central-difference solution. The open symbols include a wide variety of 
methods taken from reference 9. They include upwind schemes, central schemes, and Godunov schemes. 
The current calculations are shown as filled symbols with the eighth-order scheme flagged. The variable 
quantity is the diffusion number which changes the character of the solution from central-difference type 
to Godunov type to classical upwind type. For each value of the Cn, the optimal Dn is different. This 
optimum is much smaller for the eighth-order method. This figure shows that with a suitably chosen 
Dn, based on the given viscosity coefficient, highly accurate nonoscillatory simulations may be 
obtained. 


CONCLUSIONS 


A new method for generating difference equations approximating partial differential equations has 
been developed. The method is based on the use of a DDR approximation to the amplification factor 
that transforms the amplitude and phase information embedded in the differential operator to a finite- 
difference algorithm. While such difference equations could be derived using conventional methods, the 
use of the DDR simplified the process considerably. The difference equations are of the same order of 
accuracy as the dispersion relation associated with the differential operator, and are of high local order 
only if the coefficients are constant. 

The fourth-order method is exact for Cn of 1 and 2 and the eighth-order method is exact for Cn of 
1, 2, 3, and 4. The convection algorithms possess no dissipation and the only remaining issue is the 
phase shift due to dispersion. For problems dominated by phase accuracy, values of Cn greater than 
about 2 (for the fourth-order method) or 4 (for the eighth-order method) cause unacceptable phase dis- 
tortion. More complex algorithms based on these concepts should be useful for such problems as acou- 
stic wave diffraction, long range wave propagation through the atmosphere, the evolution of weakly 
nonlinear waves, or studies of flow stability caused by propagating waves. 
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Appendix 


This Appendix uses the symbolic algebra program Mathematica to derive 
the coefficients for the discrete dispersion relation (DDR) corresponding to 
fourth and eighth order convection. 

1. Fourth order convection. 

The following statements set up the series solution after initiating files and 
declaring variables. 

«ReIm.m 

{Re, Ira) 

«Trigonomet ry . m 

kh /: RealQ [kh]=True;bO /: RealQ [bO] =True; 
bl /: RealQ [bl]=True; b2 /: RealQ [b2] =True; 

Set up the expression for the amplification factor in eq (9) and its argument 
in eq (10) 

betap = (bO+bl Exp[-I kh] +b2 Exp[ I kh] ) / 

(b0+b2 Exp [ -I kh] +bl Exp[ I kh] ) 

-I kh I kh 

bO + E bl + E b2 


I kh -I kh 

bO + E bl + E b2 

A=Re [ComplexToTrig [Numerator [betap] ]]; 

B=Im[ComplexToTrig [Numerator [betap] ] ] ; 

Tanarg = B/A 

bl Sin [-kh] + b2 Sin [kh] 
bO + bl Cos [-kh] + b2 Cos[kh] 

Generate a Maclaurin series to fourth order with the constraint shown. Note 
that only odd powers appear. 
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sl=Series [Tanarg, {kh,0,4}] / ,bO+bl+b2->l 

bl b2 -bl b2 3 5 

(-bl + b2) kh + ( (-bl + b2) ( )) kh + 0[kh] 

6 6 2 2 

Generate the Maclaurin series for the exact dispersion relation in eq(ll). 

s2=Series [Tan[- Cn kh/2] , {kh, 0, 4} ] 

3 3 

-(Cn kh) Cn kh 5 

+ o [kh] 

2 24 

Match the coefficients of the first and third powers and solve for bl and b2. 

Coefficient [Normal [si] , kh] - 
Coefficient [Normal [s2] , kh] =0 

Cn - 

bl + b2 -- 0 

2 

Coefficient [Normal [si] , kh A 3] - 
Coefficient [Normal [s2] , kh A 3] =0 

3 

Cn bl b2 -bl b2 

+ (-bl + b2 ) ( ) »• 0 

24 6 6 22 

Solve [{%/%%) , {bl / b2}] 

2 2 
2 + 3 Cn + Cn 2 - 3 Cn + Cn 

{{bl -> , b2 -> }} 

12 12 


Clean up the expressions for bl; b2 and solve for bO from the constraint 

t=%[[l]] 

. 2 2 
2 + 3 Cn + Cn 2 - 3 Cn + Cn 

{bl -> , b2 -> } 

12 12 

bl /. t [ [1] ] ; 
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bl=% 


2 

2 + 3 Cn + Cn 
12 

b2 /. t [ [ 2 ] ] ; 
b2=% 

2 

2 - 3 Cn + Cn 
12 

bO=Factor [Simplify [l-bl-b2] ] 

-( (-2 + Cn) (2 + Cn) ) 

6 

bl=Factor[bl] 

(1 + Cn) (2 + Cn) 

12 

b2=Factor [b2] 

(-2 + Cn) (-1 + Cn) 

12 

Use these values of bO, bl, and b2 to form the 3-space-point/2-time-point 
stencil in eq (12) 

2. Eighth order convection. 

Follow exactly the same procedure as above, but include more points (and 
unknowns) in the DDR. 


b0=. ;bl=. ;b2=. ; 

kh /: RealQ[kh]=True;bO /: RealQ[bO]=True 
bl /: RealQ[bl]=True; b2 /: RealQ[b2]=True; 
b3 /: RealQ[b3]=True; b4 /: RealQ[b4]=True; 
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bet ap= (bO+bl *Exp [ - 1 kh]+b2*Exp[I kb] +b3*Exp [-2 I kh] 
+b4*Exp[2 I kh])/(bO+b2*Exp[-I kb] +bl*Exp[I kh] 
+b4*Exp[-2 I kh]+b3*Exp[2 I kh] ) 

-I kh I kh -2 I kh 2 I kh 

bO + E bl + E b2 + E b3 + E b4 


I kh -I kh 2 I kh -2 I kh 

bO + E bl + E b2 + E b3 + E b4 

A=Re [ComplexToTrig [Numerator [betap] ] ] ; 

B=Im [ComplexToTrig [Numerator [bet apj ] ] ; 
tanarg= B/A 

b3 Sin[-2 kh] + bl Sin[-kh] + b2 Sin[kh] + b4 Sin[2 kh] 

bO + b3 Cos [-2 kh] + bl Cos[-kh] + b2 Cos[kh] + b4 Cos [2 kh] 

Generate a Maclaurin series to eighth order with the constraint shown. The 
complete expression is suppressed for clarity. Only odd powers of kh will 
appear. 

sl=Series [tanarg, [kh, 0, 8 } ] /. bQ+bl+b2+b3+b4->l; 

Generate the eighth order Maclaurin series from eq (11). 

s2=Serie»[Tan[- Cn kh/2] , {kh, 0, 8} ] 

3 3 5 5 7 7 

-(Cn kh) Cn kh Cn kh 17 Cn kh 9 

+ o [kh] 

2 24 240 40320 

Match the coefficients of the first, third, fifth, and seventh powers and solve 
for bl, b2, b3, and b4 

Coefficient [Normal [si] , kh] - 
Coefficient [Normal [s2] , kh] =0 

Cn 

bl + b2 - 2 b3 + 2 b4 == 0 

2 
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Coefficient [Normal [si] , kh A 3] - 
Coefficient [Normal [s2] , kh A 3] — 0 

3 

Cn bl b2 4 b3 4 b4 
24 6 6 3 3 

-bl b2 

(-bl + b2 - 2 b3 + 2 b4) ( 2 b3 - 2 b4) == 0 

2 2 

Coefficient [Normal [si] , kh A 5] - 
Coefficient [Normal [s2] , kh A 5]=0 

5 

Cn bl b2 4 b3 4 b4 

240 120 120 15 15 

-bl b2 2 bl b2 2 b3 2 b4 

(( 2 b3 - 2 b4) - (— + — + + )) 

22 24 24 3 3 

(-bl + b2 - 2 b3 + 2 b4) - 

-bl b2 bl b2 4 b3 4 b4 

( 2 b3 - 2 b4) ( + ) == 0 

2 2 6 6 3 3 
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Coefficient [Normal [si] , kh A 7] - 
Coefficient [Normal [s2] , kh A 7]=0 

7 


17 Cn 

bl 

b2 

8 b3 

8 b4 

40320 

5040 

5040 

315 

315 


-bl 

b2 

4 b3 

4 b4 


(-(-< — 








720 

720 

45 

45 


-bl b2 

( 2 b3 - 2 b4 ) 

2 2 

bl b2 2 b3 2 b4 

( __ + — + + )) + 

24 24 3 3 

-bl b2 2 

( ( 2 b3 - 2 b4) - 

2 2 

bl b2 2 b3 2 b4 

{ __ + — + + )) 

24 24 3 3 

-bl b2 

( 2 b3 - 2 b4) ) 

2 2 


(-bl + b2 - 2 b3 + 2 b4) - 


-bl b2 


-bl 

b2 

4 b3 4 b4 

( — 

- 2 b3 

- 2 b4 ) ( — 

+ 

+ ) + 

2 2 


120 

120 

15 15 

-bl 

b2 

2 

bl 

b2 2 b3 2 b4 

( ( 

— - 2 

b3 - 2 b4 ) - 

(— + 

— + + 

2 

2 


24 

24 3 3 


bl b2 4 b3 4 b4 
6 6 3 3 
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Solve [{%, %%, %%%, %%%%}, {bl,b2,b3,b4}J 

2 3 

96 + 80 Cn + 10 Cn - 5 Cn 


{{bl -> 


420 


4 

Cn 



2 3 4 

80 Cn + 10 Cn + 5 Cn - Cn 

f 

420 


2 3 4 

24 + 50 Cn + 35 Cn + 10 Cn + Cn 

b3 -> , 

1680 


2 3 4 

24 - 50 Cn + 35 Cn - 10 Cn + Cn 

b4 -> } } 

1680 


Clean up the expressions for bl to b4 and solve for bO from the constraint. 

t=% [[ill; 
bl /. tt[l]]; 
bl=% 


2 3 4 

96 + 80 Cn + 10 Cn - 5 Cn - Cn 


420 


b2 /. t [ [ 2 ] ] ; 
b2=% 

2 3 4 

96 - 80 Cn + 10 Cn + 5 Cn - Cn 
420 

b3 /. t [ [3] ] ; 
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b3=% 


2 3 4 

24 + 50 Cn + 35 Cn + 10 Cn + Cn 
1680 

b4 /. t [ [ 4 ] ] ; 
b4=% 

2 3 4 

24 - 50 Cn + 35 Cn - 10 Cn + Cn 
1680 

bO=Factor [ Simplify [ 1 -bl -b2 -b3-b4 ] ] 

(-4 + Cn) (-3 + Cn) (3 + Cn) (4 + Cn) 

280 

bl=Factor [bl ] 

-((-4 + Cn) (2 + Cn) (3 + Cn) (4 + Cn) ) 

420 

b2=Factor[b2] 

-((-4 + Cn) (-3 + Cn) (-2 + Cn) (4 + Cn) ) 

420 

b3=Factor [b3] 

(1 + Cn) (2 + Cn) (3 + Cn) (4 + Cn) 

1680 

b4=Factor[b4] 

(-4 + Cn) (-3 + Cn) (-2 + Cn) (-1 + Cn) 

1680 

/ 

Use these values of bO - b4 to form the 5-space-point/2-time-level stencil 
shown in eq (13). This is a direct extension of the 
3-space-point/2-time-level formula shown in eq (7). 
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(LAX-WENDROFF) 


3rd ORDER UPWIND 
(LEONARD) 




Figure 1- Amplification factors for three algorithms used to compute wave propagation. 
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C N = 0.3 




Figure 3 - Phase accuracy for dissipation-free second-, fourth-, and eighth-order algorithms. 
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Figure 4 — Propagation of a pulse in an infinite medium. Cn = 1.23. (a) Initial (A) and final positions 
with theoretical curves superimposed (B); (b) amplitude spectra of initial and final pulses; (c) phase 
dispersion for eighth-order scheme; (d) phase dispersion for fourth-order scheme. Exact phase shown 
dashed. 
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Figure 5 - Propagation of a pulse in an infinite medium. Cn = 4.78. (a) Initial (A) and final positions 
with theoretical curves superimposed (B); (b) amplitude spectra of initial and final pulses; (c) phase 
dispersion for eighth-order scheme; (d) phase dispersion for fourth-order scheme. Exact phase shown 
dashed. 
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Figure 8.- Amplitude spectra and phase characterisitics for the propagation of a sharp pulse with 

viscosity. Cn =4 .78, Dn = 0.00475. 
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Figure 9.— Phase accuracy of the discrete dispersion relation. Cn = 0.5. (a) Amplitude spectrum for a 
step function; (b) eighth-order scheme; (c) fourth-order scheme; (d) second-order scheme. 
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Figure 10.- Propagation of a step after 45h space intervals. Cn - 0.5. (a) Second-order scheme; 
(b) fourth-order scheme; (c) eighth-order scheme. Initial and Final position shown dashed. 
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Figure 11- Effect of diffusion number on step function propagation using fourth-order algorithm. 

Cn = 0.5. Initial and final position shown dashed. 
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Figure 13 - Effect of diffusion number on the amplitude spectrum for the propagation of a step function. 
Cn = 0.5. Eighth-order algorithm. Solid line: spectrum at initial instant. Dashed line: spectrum after 
45h intervals with no viscosity. Dotted line: computed spectrum after 45h intervals. 
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Figure 14.- Accuracy of step function calculations. Open symbols: Leonard (ref. 9). Closed symbols: 

current simulations. 
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